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Outline 


An implicit, two-factor, split flux, finite volume Euler equations solution 
algorithm is applied to the time-accurate solution of transonic flow about an 
NACA0012 airfoil and a rectangular planform supercritical wing undergoing pitch 
oscillations. Accuracy for Courant numbers greater than one is analyzed. Freezing 
the flux Jacobians can result in significant savings for steady-state solutions; the 
accuracy of freezing flux Jacobians for unsteady results is investigated. The Euler 
algorithm results are compared with experimental results for an NACA 0012 and a 
rectangular planform supercritical wing (Figure 1). 


• ALGORITHM 

• TIME ACCURACY FOR COURANT NUMBERS GREATER THAN ONE 

• FREEZING FLUX JACOBIANS 

• COMPARISON WITH EXPERIMENT 

• NACA0012 

• RECTANGULAR SUPERCRITICAL WING 

Figure 1 
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Three Dimensional Euler Equations 


The Euler equations model inviscid flow with entropy and vorticity being gener- 
ated by shocks. Many unsteady flows of practical interest require that viscous 
effects be included. One advantage of the Euler equations solution algorithm is that 
it is easily extended to include viscous effects. 

The three-dimensional time-dependent Euler equations in conservation form are 
written in a general time-dependent boundary conforming curvilinear coordinate 
system. This time-dependent coordinate transformation provides for a dynamic grid 
that can follow the motion of the body (Figure 2). The details of this transforma- 
tion are given in Reference 1. 
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Figure 2 
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Algorithm 


An implicit split flux finite volume solution algorithm for the Euler equations 
is obtained. Flux Jacobians with superscript L correspond to information 
propagating in the positive curvilinear coordinate direction, and flux Jacobians 
with superscript R correspond to information propagating in the negative curvilinear 
coordinate direction. The equation is then factored into two factors as discussed 
by Whitfield in Reference 2. The first equation of the two-pass scheme requires the 
solution of a sparse lower block triangular system by a simple forward 
substitution, while the second equation requires the solution of a sparse upper 
block triangular system by a simple backward substitution. 

This algorithm is first-order accurate in time and second order accurate in 
space. A simple modification to use three point backward time differencing will 
result in second order time-accuracy (Figure 3), 


• FINITE VOLUME 

• FLUX SPLIT 

• IMPLICIT 
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WHERE 

R n = 8f F n 4* 8-q e G n + <5£ e H n 
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Figure 3 
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ORIGINAL PAGE IS 

NACA0012 221 x 20 'C' Grid 0F P00R QUALITY 


All of the Euler computations for the NACA0012 airfoil used the 221 x 20 r C f 
algebraic grid shown below. The grid was generated using Joe Thompsons grid 
generation code (Reference 3). See Figure 4. 



NACA0012 221 x 20 'C Grid 


Figure 4 
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NACA0012 Lift Coefficients 


ORIGINAL PAGE IS 
OF POOR QUALITY 


Calculations were made for a NACA0012 airfoil at Mach = 0.755 oscillating in 
pitch about the 1/4 chord point with a reduced frequency of k = (<^c/V ) = 0.1628. 

The airfoil had a mean angle of attack = 0.016 degrees and an unsteady alpha ■ 2.51 
degrees. Steady state solutions were obtained at the mean conditions prior to an 
abrupt start of the oscillatory motion. The calculations were performed for four 
complete cycles of motion. The figure below shows lift coefficient vs time for 
three different time step sizes. DTMIN - 0.01 gave 5000 time steps per cycle of 
motion and corresponds to a maximum Courant number of 10. DTMIN = 0.10 was 500 time 
steps per cycle of motion and corresponds to a maximum Courant number of 100. DTMIN 
= 0.20 was 250 time steps per cycle of motion and corresponds to a maximum Courant 
number of 200. Lift coefficients were only slightly different for the various time 
step sizes (Figure 5). 


11=0.755, k=0 . 1 628 , MEAN ALPHA=0.Q16, UNSTEADY ALPHA=2.51 
221X20 ALCEBRAIC ‘C’ GRID 
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NACA0012 Unsteady Pressure Distributions 
340° of Oscillatory Motion 


After 340° of oscillatory motion, the angle of attack is increasing through 
-0.84 degrees. A shock has formed on the lower surface near the 44% chord location. 
At this point in the motion, the shock is nearly stationary and of maximum strength. 
The figure below shows the coefficient of pressure distribution along the airfoil. 
Calculations were done using both first and second order time-accurate differencing 
and various time step sizes. The results show that all the methods are nearly 
equivalent for this case (Figure 6). 


M=0.755, Unsteady Alpha=2.51, Mean Alpha=0.016, k=0.l628, 221X20 'C # Grid 

Two Pass Algorithm 



340° of Oscillatory Motion, Increasing Through 0-84° 


Figure 6 
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NACA0012 Unsteady Pressure Distributions 
25° of Oscillatory Motion 


Shortly after the next cycle of motion begins, the shock on the lower surface 
starts to collapse and move forward. As the angle of attack increases through 
1.09°, the shock speed becomes maximum. The figure below shows coefficient of 
pressure along the airfoil for first and second order differencing and the same time 
step sizes as in Figure 5. Very little difference is observed between first and 
second order solutions, but considerable difference in the shock location is noted 
for different time step sizes. However, as noted previously, the difference in lift 
coefficient is relatively small (Figure 7). 


M-0 . 755 , Unsteady Alpha=2.5I. Mein Alpha=0.016. k=0.1628, 221X20 *C' Grid 

Tvo Pass Algorithm 



25 ° of Oscillatory Motion, Increasing Through 1 - 09 ° 


Figure 7 
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Shock Location 


Defining the shock location as the axial location where the pressure coefficient 
equals the critical pressure (for Mach * 0.755, C * ■ -0.5143), the figure below 
shows shock location given at different time step F sizes. The figure shows that the 
shock location appears asymptotic to a value as the time step decreases. This 
indicates that reduction of the time step size below 0.01 should not be expected to 
change the shock location appreciably (Figure 8). 
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Freezing the Flux Jacobians for Steady State Conditions 


The implicit equations to be solved have the flux Jacobians as coefficients* 
Since the flux Jacobians, A^, A R , B^, etc., are functions of the current values of 
Q, they should be updated at each time step. However, D. L. Whitfield has shown 
steady calculations in which he did not update (froze) the Jacobians and at 
convergence obtained identical results with calculations updating the Jacobians each 
time step. The table below shows the obvious computational savings by not doing the 
extra calculations each time step (Figure 9). 
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Figure 9 
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Freezing the Flux Jacobians for Unsteady Calculations 


The effect of freezing the Jacobians for unsteady calculations were studied 
using the oscillatory NACA0012 as the test case. The conditions tested were the 
same as previously shown. The table below again shows the obvious savings from 
updating every 10th step (which equates to 3.6° of oscillatory motion) and for 
never updating the Jacobians. The case listed as never updated used Jacobians from 
the steady state condition just prior to start of motion. Other cases tried were 
Jacobians from the freestream starting conditions, updating every 23th step and 
updating every 50th step. Each of these resulted in stability problems (Figure 10). 


JACOBIAN UPDATE 
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NACA0012 PITCHING ABOUT 1/4 CHORD 


Figure 10 
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Effect of Frozen Jacobians on Lift Coefficient 


Lift coefficient vs time for each of the two frozen Jacobian cases which proved 
to be stable are compared with the solution from updating every time step. All 
three solutions used 1st order time-accurate algorithm and a time step size of 0.05. 
One cannot discern a difference in the three curves shown below (Figure 11). 


M=0 . 755 , Unsteady Alpha=2.51, Mean Alpha=0 . 016 , k=0 . 1628, DT=0 . 05 



TIME (DEG) 


Figure 11 
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NACA0012 Unsteady Pressure Distributions 
25 Degrees 


A plot of pressure coefficient along the body for the three cases at 25 degrees 
of oscillatory motion is, shown below in Figure 12. This is the same flow condition 
described earlier, when the angle of attack is increasing through 1.09 degrees and 
the shock speed is near a maximum as the shock on the lower surface collapses. Very 
little difference is noted in the three curves with the only perceivable difference 
being in the shock location similar to the results shown in Figure 7. 



Figure 12 
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NACA0012 Unsteady Pressure Distributions 
70 Degrees 


A plot of pressure coefficient along the body for the three cases at 70 degrees 
of oscillatory motion is shown below. This condition is when angle of attack is 
increasing through 2.37 degrees and the shock on the upper surface is strengthening 
and moving downstream. Again, the three solutions are identical except at the shock 
location. Note that the never updated Jacobian shows a slight ringing action near 
the shock (Figure 13). 


M=0 . 755 , Unsteady Alpha=2.51, Mean A 1 pha=0 .016, k=0 . 1628 , 22 1X20 ' C ' Grid 
Two Pass Algorithm, DT=0.05, lst-Order in Time 



Figure 13 
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NACA0012 Lift Coefficients 


A comparison of calculated lift coefficients with the experimental results of 
Landon (Reference 4) is given in Figure 14. The experimental results were given for 
one cycle of motion and are duplicated through four cycles on the figure to compare 
to the calculations. Examination of experimental lift coefficient reveals a bias 
towards positive lift not consistent with the small mean angle of attack reported 
for the symmetric airfoil. A correction to the mean angle of attack was calculated 
to account for this bias. The following figures show calculations using both the 
nominal angle of attack of 0.016 degrees and th£ 'corrected' angle of attack of 
0.375 degrees. The freestream Mach is 0.755, reduced frequency is 0.1628, and the 
unsteady angle of attack amplitude is 2.51 degrees. 


M=0 . 755 , k=0 . 1628, UNSTEADY ALPHAS. 51 
221120 ALGEBRAIC 'C' GRID 



Figure 14 
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NACA0012 Moment Coefficients 


Moment coefficients do not compare with experiment as well as lift, particularly 
for the negative moments. The 'corrected' mean alpha of 0.375 did not improve the 
comparison (Figure 15). 


M=0 . 755 , k=0 . 1628, MEAN ALPHA=0.016, UNSTEADY ALPHA=2.51 
221X20 ALGEBRAIC 'C' GRID 



Figure 15 
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ORIGINAL PAGE IS 
OF POOR QUALITY 

NACA0012 Unsteady Pressure Distributions 


The comparison of calculated and experimental pressure distributions is shown in 
Figures 16 through 21. 


M=0 . 755 , Unsteady Alpha=2.5l, Mean Aipha=0.016, k=O.l620, 221X20 'C' Grid 



M=0 . 755 , Unsteady Alpha=2.51, Mean Alpha=0.016, k=0 . 1620 , 221X20 'C* Grid 



Figure 17 
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NACA0012 UNSTEADY PRESSURE DISTRIBUTIONS 

M=0 . 755 , Unsteady Alpha=2.51, Mean Alpha=0 . 0 16 , k=0.1628. 221X20 'C' Grid 



NACAOO 1 2 UNSTEADY PRESSURE DISTRIBUTIONS 

M=0 . 755 , Unsteady Alpha=2.51, Mean Alpha=0 . 0 16 , k=0.1628, 221X20 'C' Grid 



Figure 19 





ORIGINAL PAGE IS 
OF POOR QUALITY 


NACA0012 UNSTEADY PRESSURE DISTRIBUTIONS 

M=0 . 755 , Unsteady Alpha=2.51, Mean Alpha=0.018, k=0.1628, 221X20 'C' Grid 



NACAOO 1 2 UNSTEADY PRESSURE DISTRIBUTIONS 

11=0.755, Unsteady Alpha=2.51, Mean Alpha=0.016, k=0.1628, 221X20 ’C' Grid 



Figure 21 




Rectangular Supercritical Wing 


Euler calculations were performed to compare to the experimental results of 
Ricketts, et al (Reference 5). Calculations at a higher reduced frequency, k = 
0*714, compare similarly to those presented here, k = 0.358. At higher Mach 
numbers, however, the comparison was much poorer due to mis location of the shock by 
the inviscid Euler code. The time step size used resulted in 360 time steps per 
cycle of motion for the k = 0.358 case. Maximum Courant numbers near 500 occurred 
in the vicinity of the wing tip for this time step size. Fourier analysis of the 
third cycle of oscillation yielded the magnitude and phase of the unsteady pressures 
shown on subsequent figures. The three cycles of motion used 6357 seconds on a CRAY 
X-MP (Figure 22). 


• EXPERIMENT BY RICKETTS* SANDFORD, WATSON* AND SEIDEL 
NASA TM 85765* AUG 84 

• RECTANGULAR PLANF0RM* ASPECT RATIO 4 (FULL SPAN) 

• 121 THICK SUPERCRITICAL AIRFOIL 

• OSCILLATORY PITCH ABOUT 46Z CHORD 

• CONDITIONS 

• MACH 0.70 

• 4 DEGREES fCAN ALPHA 

• ONE DEGREE UNSTEADY ALPHA 

• REDUCED FREQUENCY = 0-358 = 


Figure 22 


234 



ORIGINAL PAGE IS 

Rectangular Wing Blocked Grid OF POOR QUALITY 


The wing calculations were carried out on a grid broken into four blocks to 
obtain the solution using only 2.6 million words of memory. The entire grid has 
dimensions 101x25x27. Block I contains all points below the wing. Block II contains 
points wrapping around the wing tip. Block III contains all points above the wing, 
and Block IV contains all points downstream of the wing (Figure 23). The method used 
to obtain time-accurate solutions on blocked grids is described in Reference 6. 

Motion of the wing is modelled by pitching the entire grid containing the wing as a 
rigid body using the time-dependent coordinate transformation described earlier. 
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Rectangular Wing Pressures 
at 60% Semispan 


Steady pressures compare fairly well. There is a slight overshoot at the 
leading edge and the captured shock is aft of the experimental location (Figure 24). 
Unsteady pressure magnitude compares well except that the shock spike is downstream 
of the experimental location. Unsteady pressure phase calculations show excellent 
agreement up to the highly cambered trailing edge where the experiment and 
calculation differ slightly (Figures 25 and 26). 


RECTANGULAR SUPERCRITICAL WING 

STEADY PRESSURE DISTRIBUTIONS 
Mach 0.7.01; 4.0 degrees Angle of Attack 



Figure 24 
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RECTANGULAR SUPERCRITICAL WING 

UNSTEADY PRESSURE DISTRIBUTIONS 
Mach 0.7; 4 deg Steady Alpha; 1 deg Usleady Alpha 



Magnitude at 60% Semi span 

Figure 25 


RECTANGULAR SUPERCRITICAL W 

UNSTEADY PRESSURE DISTRIBUTIONS 
Mach 0.7; * deg Steady Alpha; 1 deg Unsteady 
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Phase at 60% Semspan 


Figure 26 
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Rectangular Wing Pressures 
at 95% Semispan 


Steady pressures shown in Figure 27 underexpand slightly in the leading edge 
region but follow the experimental results very closely thereafter. The peak in 
unsteady pressure magnitude shown in Figure 28 near the leading edge is likewise 
calculated to be smaller than experiment. Phase results shown in Figure 29 are in 
excellent agreement. 


RECTANGULAR SUPERCRITICAL WING 

STEADY PRESSURE DISTRIBUTIONS 



Figure 27 
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PHASE MAG 


RECTANGULAR SUPERCRITICAL WING 

UNSTEADY PRESSURE DISTRIBUTIONS 
Mach 0.7; 4 de( Steady Alpha; l deg Usteady Alpha 



Figure 28 


RECTANGULAR SUPERCRITICAL WING 

UNSTEADY PRESSURE DISTRIBUTIONS 
Mach 0.7; 4 deg Steady Alpha; 1 deg Unsteady Alpha 



Phase at 95% Semispan 


Figure 29 





CONCLUSIONS 


• TIME-ACCURACY CAN BE OBTAINED AT 
COURANT NUMBERS MUCH GREATER THAN ONE 


• FOR THE CASE PRESENTED, FREEZING FLUX 
JACOB IANS HAD LITTLE EFFECT ON TIME-ACCURACY 


• EULER CALCULATIONS COMPARE WELL WITH NACA0012 
AND SUPERCRITICAL WING EXPERIMENT 


• VISCOSITY REQUIRED TO ACCURATELY MODEL 
SUPERCRITICAL WING 


Figure 30 
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